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STABILITY OF DISTANT SATELLITES OF THE GIANT PLANETS IN THE SOLAR SYSTEM 
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ABSTRACT 

We conduct a systematic survey of the regions in which distant satellites can orbit stably around the four giant 
planets in the solar system, using orbital integrations of up to 10 9 yr. In contrast to previous investigations, we 
use a grid of initial conditions on a surface of section to explore phase space uniformly inside and outside the 
planet's Hill sphere (radius th; satellites outside the Hill sphere sometimes are also known as quasi- satellites). 
Our confirmations and extensions of old results and new findings include the following: (i) many prograde and 
retrograde satellites can survive out to radii ~ 0.5rn and ~ 0.7rn, respectively, while some coplanar retrograde 
satellites of Jupiter and Neptune can survive out to ~ rn; (ii) stable orbits do not exist within the Hill sphere 
at high ecliptic inclinations when the semi-major axis is large enough that the solar tide is the dominant non- 
Keplerian perturbation; (iii) there is a gap between ~ rn and 2rn in which no stable orbits exist; (iv) at distances 
> 2rn stable satellite orbits exist around Jupiter, Uranus and Neptune (but not Saturn). For Uranus and Neptune, 
in particular, stable orbits are found at distances as large as ~ 10rn; (v) the differences in the stable zones 
beyond the Hill sphere arise mainly from differences in the planet/Sun mass ratio and perturbations from other 
planets; in particular, the absence of stable satellites around Saturn is mainly due to perturbations from Jupiter. 
It is therefore likely that satellites at distances > 2rn could survive for the lifetime of the solar system around 
Uranus, Neptune, and perhaps Jupiter. 

Subject headings: celestial mechanics - planets and satellites: formation - planets and satellites: general - 
minor planets, asteroids 



1. INTRODUCTION 

Most of the satellites of the four giant planets in the solar 
system can be divided into two groups, usually called the reg- 
ular and irregular satellites. Regular satellites orbit close to 
the planet (within ~ 0.05rn, where rn is the Hill radius 1 ), and 
move on nearly circular, prograde orbits that lie close to the 
planetary equator. Irregular satellites are found at distances 
~ 0.05rH-0.6rH, with large orbital eccentricities and incli- 
nations, on both prograde and retrograde orbits. An alterna- 
tive division between regular and irregular satellites is given 
by the critical semi-major axis (e.g., Goldreich 1966; Burns 
1986), a CY i t = (2fiJ2R 2 a 3 p ) 1 / 5 ; those with a > a cr i t are classified 
as irregular satellites. Here J 2 is the planet's second zonal har- 
monic coefficient (augmented by any contribution from the in- 
ner regular satellites) and R is the planet's radius. This critical 
radius marks the location where the precession of the satel- 
lite's orbital plane is dominated by the Sun rather than by the 
planet's oblateness. The current number ratios of irregular to 
regular satellites are 55/8 for Jupiter, 35/21 for Saturn, 9/18 
for Uranus, and 7/6 for Neptune (e.g., Jewitt & Haghighipour 
2007). The regular satellites are likely to have formed within 
a circumplanetary disk of gas and solid bodies. The kine- 
matic differences between regular and irregular satellites sug- 
gest that the latter must have formed through a quite different 
mechanism, most likely capture from the circumstellar disk 
(for a recent review, see Jewitt & Haghighipour 2007). 

The search for irregular satellites of the giant planets has 
been fruitful in recent years, owing mainly to modern high- 
sensitivity, large-scale CCDs (e.g., Gladman et al. 1998, 2000, 
2001; Holman et al. 2004; Kavelaars et al. 2004; Sheppard 
& Jewitt 2003; Sheppard et al. 2005, 2006). Up to 2007, 
106 irregular satellites of the giant planets had been discov- 

1 The Hill radius is defined as rn = ^(/i/3) 1 / 3 , where a p is the semi-major 
axis of the planet orbit and ji = m p /(m p +M@) with m p the planet mass. 



ered, compared to 53 regular satellites. Two features stand 
out in the distributions of orbital parameters of these irreg- 
ular satellites. First, retrograde irregular satellites extend to 
larger semi-major axes than prograde ones (~ 0.6rn compared 
to ~ 0.4rn); second, satellites with orbital inclination in the 
range ~ 60° - 130° relative to the ecliptic are absent. 

A number of authors have shown that these features can be 
explained reasonably well by the requirement that the satellite 
orbits be stable. Henon (1969, 1970) studied the planar circu- 
lar restricted three-body problem in Hill's (1886) approxima- 
tion, where the mass ratio /i — > while the radii of interest 
shrink to zero as /i 1 / 3 . He showed that prograde satellite or- 
bits are stable up to a mean distance from the planet ~ 0.4r H , 
while retrograde satellite orbits can be stable at much larger 
distances from the planet. Thus it is not surprising that ret- 
rograde satellites are found at larger distances than prograde 
ones. Hamilton & Krivov (1997) studied the dynamics of 
distant satellites of asteroids in heliocentric orbits using a 
"generalized Tisserand constant" and, among other conclu- 
sions, confirmed that retrograde orbits are more stable than 
prograde ones. Carruba et al. (2002) used a combination of 
analytic arguments and numerical integrations to show that 
high-inclination orbits inside the Hill sphere exhibit large ec- 
centricity oscillations (Kozai oscillations; Kozai 1962) due to 
secular solar perturbations. They found that orbits with incli- 
nations (relative to the planetary orbital plane) between 55° 
and 130° are generally unstable, thus explaining the absence 
of irregular satellites on high-inclination orbits. Nesvorny et 
al. (2003) performed detailed orbital integrations of the four 
giant planets plus a grid of test-particle satellites for intervals 
of 10 6 -10 8 yr. They confirmed that retrograde satellites can 
be stable at larger radii than prograde ones, and that highly in- 
clined orbits are unstable. They argued that the largest semi- 
major axes at which satellites of the four giant planets could 
survive for times comparable to the lifetime of the solar sys- 
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tern were ~ 0.7r H for retrograde satellites and ~ 0.4r H for 
prograde ones, and that these upper limits were achieved only 
for nearly circular orbits close to the plane of the ecliptic. 

Other authors have examined the possibility that stable 
satellite orbits exist with mean distance from the planet > rn. 
In the planetocentric frame, the dominant force on such satel- 
lites is due to the Sun, rather than the planet 2 . Neverthe- 
less, the satellite remains close to the planet because it is in 
a 1:1 resonance in the sense that its heliocentric mean lon- 
gitude librates around that of the planet; the resulting orbit 
relative to the planet is a retrograde ellipse with axis ratio 
2:1, the short axis pointing towards the Sun, and synodic pe- 
riod equal to the planet's orbital period. The analytical the- 
ory of such orbits is described by Jackson (1913), Lidov & 
Vashkov'yak (1994a,b), Mikkola & Innanen (1997), Namouni 
(1999), Mikkola et al. (2006), and others. Henon's (1970) nu- 
merical analysis of Hill's approximation to the planar circular 
restricted three-body problem suggests that stable retrograde 
satellites can exist at arbitrarily large distance from the planet. 
Benest (1971) confirmed that stable retrograde orbits at large 
distances persist in the elliptic restricted three-body problem, 
where the mass ratio and eccentricity were chosen to match 
those of Jupiter. Wiegert et al. (2000) demonstrated that ret- 
rograde satellites of Uranus and Neptune could be stable for 
up to 10 9 yr at distances up to ~ 10rn, suggesting that primor- 
dial objects of this type could still exist in the solar system 
although none are currently known. 

Despite the number and quality of these investigations, 
there are several unanswered questions that lead us to revisit 
the problem of orbital stability of satellites at large distances 
from the host planet, (i) Wiegert et al. (2000) found stable 
satellite orbits beyond the Hill radius only for Uranus and 
Neptune, not Jupiter or Saturn. What is the reason for this 
difference? The possibilities include differences in the plan- 
etary masses and orbital eccentricities, or different perturba- 
tions from neighboring planets, (ii) Wiegert et al. (2000) ex- 
plored orbits outside the Hill radius, while Nesvorny et al. 
(2003) explored orbits inside the Hill radius (indeed, in the 
former paper the integrations were terminated when the par- 
ticles entered the Hill sphere of radius rn around the planet, 
while in the latter paper the integrations were terminated when 
the particles exited the Hill sphere). Are there stable satellite 
orbits that cross the Hill sphere? (iii) As we shall describe 
further in ^2) the grids of initial conditions used by Nesvorny 
et al. and Wiegert et al. do not provide a complete exploration 
of the phase space in which stable satellite orbits exist. 

The primary goal of this paper is to map out the entire sta- 
bility region in phase space — both inside and outside the Hill 
sphere — in which satellite orbits that can survive around the 
four giant planets for times comparable to the age of the so- 
lar system (our main integrations last for up to 100 Myr). We 
describe our setup in ^2) and present the results in §|3] We 
conclude and discuss our results in §@] 

Following Fabrycky (2008), we shall define a "satellite" of 
a planet to be a small body whose distance from the planet 
never exceeds the semi-major axis of the planet, a p . This def- 
inition excludes bodies on Trojan orbits around the triangu- 
lar Lagrange points, bodies on horseshoe orbits, and objects 
such as asteroid 2003 YN107 (Connors et al. 2004), which 
oscillates between a horseshoe orbit and an orbit centered on 
Earth. This definition seems simple and reasonable to us, but 

2 Hence these are sometimes called "quasi-satellites" (Lidov & 
Vashkov'yak 1994a,b; Mikkola & Innanen 1997). 



other definitions are common in the literature. Many authors 
define "satellite" to be an object that always remains within 
the Hill sphere of the planet or whose Jacobi constant con- 
strains it to remain within the last closed zero-velocity sur- 
face around the planet. Benest (1971) defines a satellite to 
be a body whose heliocentric orbital frequency is the same as 
the planet's, but whose synodic frequency around the planet 
is non-zero. Wiegert et al. use the term "quasi- satellite" for 
an object that remains outside the Hill sphere but whose he- 
liocentric longitude difference from the planet never exceeds 
120° and regularly passes through zero. However, the term 
"quasi- satellite" is confusing because it is also used for ob- 
jects such as 2003 YN107 that spend part of their time on 
horseshoe orbits and thus are only temporarily satellites in 
our sense. 

2. METHODS 

Although all of our results are based on direct numerical 
integrations of the N-body problem (Sun, one or four giant 
planets, plus a test particle orbiting one planet), we shall find 
it useful to interpret our results in terms of the coordinates and 
notation used by Henon (1970) in the exploration of satellite 
orbits in Hill's approximation. 

2.1. Hill's approximation 

When studying satellite motions near a planet (r <C rn) it 
is conventional to employ a non-rotating planetocentric coor- 
dinate system, which we denote as (xyz). However, in Hill's 
approximation to the circular restricted three-body problem, 
it is more convenient to use a rotating planetocentric coordi- 
nate system (^rjQ, where £, r] and £ are scaled coordinates 
in the rotating frame in which the planet is at the origin, the 
£ axis is along the direction opposite the Sun and the ( axis 
is perpendicular to the Sun-planet orbital plane. In Hill's for- 
mulation the unit of length is \± l ^a v , and the unit of time is 
n~ l where n = [G(M +m p )/a 3 p ] 1 / 2 is the mean motion of the 
planet. As usual, the orbit of the planet in the inertial frame 
is counter-clockwise as viewed from the positive z or £ axis. 
In Hill's coordinate system the collinear Lagrangian points L\ 
and L 2 are located at r\ = ( = 0, £ = ±3 -1 / 3 ~ 0.6934, and the 
Hill radius is rn = 3 -1 / 3 . Similar definitions are used in this 
paper when the planet orbit is eccentric and/or perturbed by 
other planets; in this case the £ axis points away from the in- 
stantaneous position of the Sun, the ( axis is perpendicular to 
the instantaneous orbital plane of the planet around the Sun, 
and a p is the initial semi-major axis of the planet. 

In the circular restricted three-body problem, Hill's approx- 
imation is achieved by taking the limit /i — > 0, where the equa- 
tions of motion reduce to (e.g., Henon 1974; Murray & Der- 
mott 1999): 

^ 2 ^ +3 ^~(£ 2 +r/+c 2 ) 3/2 ' " i= ~ 2 ^~ (ew+c 2 ) 3/2 ' 



(£ 2 + 7? 2 + C 2 ) 3/2 ' 

There exists an integral of motion for these equations, 

r = ^ + ( g 2 + ??2 2 +C2)1/ 2 -C 2 -^ 2 ^ 2 + C 2 ), (2) 

which corresponds to the Jacobi constant in the circular re- 
stricted three-body problem. 

For the moment let us restrict ourselves to motion in the 
Sun-planet orbital plane, so ( = ( = at all times. Then to 
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Fig. 1 . — Sampling of initial conditions in terms of the Henon diagram. The gray-shaded regions are forbidden. The Lagrange points L\ and L2 are marked by 
*. The horizontal dashed line at £ = separates retrograde orbits from prograde ones (as defined in the rotating planetocentric frame). Left: The region shaded 
by vertical lines is an approximate reproduction of the stable region as estimated in Henon (1970). The dotted region is where the osculating Kepler elements 
correspond to bound elliptical orbits (a > 0, e < 1). Right: Curved lines represent the initial conditions derived using osculating elements in the high-resolution 
survey of Nesvorny et al. (2003), color-coded according to eccentricity. Solid and long-dashed lines represent orbits which are prograde and retrograde in 
the non-rotating planetocentric frame respectively. There are two sets of lines for each eccentricity corresponding to argument of pericenter uj = 0° and 180° 
respectively. The short segments at the lower right below L\ are extensions of the e = 0.50 and 0.75 branches which are prograde in the non-rotating frame. Thus 
orbits that are retrograde in the rotating frame can be prograde in the non-rotating frame. The initial conditions for zero-inclination orbits sampled by Wiegert et 
al. (2000) are shown as filled circles for Uranus and open circles for Jupiter. 



study orbital motions we may use a surface of section de- 
fined by 77 = 0, 77 > 0. The trajectory in the four-dimensional 
(£, 77, £, 77) phase space is then represented by a set of points in 
the (£, O plane, and for a given value of the Jacobi constant T 
the other two phase- space coordinates can be derived from 



77 = 0, 



*7=(3£ 2 -£ 2 +^-r 



1/2 



(3) 



We define "prograde" and "retrograde" in the rotating frame 
unless otherwise noted. Thus retrograde orbits have £ < in 
this surface of section and prograde orbits have £ > 0. 

A drawback of this surface of section is that a different plot 
is needed for each value of the Jacobi constant T. To obtain a 
global view of the dynamics, we use a different surface of sec- 
tion defined by 77 = 0, f = 0, 77 = (3£ 2 + 2/|f | -T) 1 / 2 . A trajec- 
tory is represented by a point in the (T, £) plane. This surface 
of section was introduced by Henon (1969; 1970), and we 
shall call it the Henon surface of section or Henon diagram. 
The Henon diagram, like any surface of section, will not show 
orbits that do not cross it; the usefulness of the Henon diagram 
derives from the observation that most stable orbits periodi- 
cally pass close to the point 77 = £ = — for example, this oc- 
curs for nearly Keplerian orbits close to the planet when their 
line of apsides precesses past the Sun-planet line. The orbits 
not shown on the Henon diagram include those confined to 
some resonant islands, which should occupy a small fraction 
of phase space, and escape orbits, which we are not interested 
in anyway. 

Fig. Q] (left panel) is a Henon diagram modeled on Fig- 
ure 12 of Henon (1970). The Lagrange points are at (I\£) = 
(3 4 / 3 ,±3 -1 / 3 ) = (4.32675, ±0.69336). Forbidden regions, in 
which rj 1 = 3£ 2 + 2/|£| -V would be negative, are shaded in 
gray. The stable regions of phase space, as estimated by 



Henon, are denoted by vertical stripes. The diagram shows 
that retrograde satellites (£ < 0) have a larger stable region 
than prograde satellites (£ > 0), a conclusion consistent with 
the numerical studies described in §Q] Moreover, the stable 
band in this diagram that begins at (T,£) = (-1.4,-1.2) and 
stretches downward to the left shows that retrograde satellites 
can be stable at distances much larger than the Hill radius; in 
fact, this band continues to arbitrarily large negative values of 
T and £ (see Figure 13 of Henon 1970), so retrograde satellites 
can be stable at arbitrarily large distances from the planet, at 
least in Hill's approximation to the planar circular restricted 
three-body problem. 

In future discussions we divide the stable regions in Fig. Q] 
into three branches: the inner prograde branch (£ > 0), the 
inner retrograde branch (£ < and T > 0), and the outer ret- 
rograde branch (£ < and T < 0). 

A simple and rather complete way to sample initial condi- 
tions in the planar three-body problem is to use the Henon 
diagram, i.e., to sample uniformly in the plane. As de- 
scribed above, this approach is based on the assumption that 
most stable orbits periodically have their apocenter or peri- 
center on the Sun-planet line. Note that even without invok- 
ing Hill's approximation the question of which initial con- 
ditions on the Henon diagram correspond to stable orbits is 
well-posed. Accordingly, we may present our stability results 
in terms of the Henon diagrams, even though our orbit inte- 
grations do not use Hill's approximation. 

We may compare this approach to the grids of initial con- 
ditions used in other investigations of the stability of satel- 
lite orbits. The initial conditions for Nesvorny et al.'s "high- 
resolution survey" were chosen from a grid of planet-centered 
osculating Keplerian orbital elements, with semi-major axis 
a given typically by a/ru = 0.1-1, eccentricity e = 0-0.75, 
inclination i = 0°-180°, argument of pericenter uj = 0°,90°, 
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and the other elements distributed uniformly between 0° and 
360°. The right panel of Fig. [T] shows similar initial condi- 
tions on the Henon diagram (i = 0° or 180° and u = 0° or 
180°). The conversions from osculating elements to (T,0 
were done using equations (8) and (10) in Henon (1970). It 
is clear that the initial conditions sampled in Nesvorny et al. 
do not provide a complete exploration of the phase space in 
which stable satellite orbits could exist; in particular, they 
completely missed the stable region that extends beyond the 
Hill sphere (of course, such orbits are also excluded from their 
study by their artificially imposed escape criterion r > rn). In 
fact, most of the stable orbits beyond the Hill sphere have hy- 
perbolic osculating elements. In Fig. [T] we plot the boundaries 
that separate regions of hyperbolic osculating elements from 
those with elliptical osculating elements, where the latter are 
shaded by a dotted pattern. The functional forms of these 
boundaries are: 



r=2e+2 3 / 2 \t\ l/2 

T = 2e + 2 3 ^ 2 



(£ < 0) , 

(0 < f < 2 1/3 ) , 

(e<-2 1 / 3 ). 



(4) 
(5) 
(6) 



The initial conditions explored by Wiegert et al. (2000) 
were chosen from a grid of heliocentric osculating Keplerian 
elements, these being the same as the elements of the host 
planets except for the eccentricity and inclination. The ec- 
centricity was typically chosen in the range e = 0-0.5 and 
inclination in the range 0°-30°. With this procedure, zero- 
inclination orbits appear in the Henon diagram along the lo- 
cus 

r = 2/|e| -e with £ = -^" 1/3 , (7) 

where the expression for Y is evaluated using Hill's approx- 
imation. The grid sampled by Wiegert et al. for i = is also 
shown in the right panel of Fig. [T] converted from the helio- 
centric frame using Hill's units (but without Hill's approxi- 
mation); for clarity, only Jupiter and Uranus are shown. Al- 
though Wiegert et al.'s initial conditions do probe the stability 
region found by Henon beyond the Hill sphere, the coverage 
is far from complete. 

To extend our study to three-dimensional orbital motions 
we use a surface of section at r] = ( = £ = 0, rj > 0. In the 
rotating frame, we define the initial inclination angle I by 



tan/ = 



(8) 



f=0 



such that the initial 77 and ( component velocities are 

1/2 



77 = cos/(3^ + — -r 



C = sin/(3£ 2 +^-r 



1/2 



(9) 



Since 1) > 0, the inclination is restricted to the range -90° < 
7<90°. 

Therefore each point in the (T, £) plane represents a unique 
set of initial conditions for a given inclination. The usefulness 
of the Henon diagram in this case is based on the assumption 
that most stable orbits periodically have their line of apsides 
and their line of nodes simultaneously on the Sun-planet line. 
This assumption is not always valid: it requires that the argu- 
ment of pericenter oj is periodically or tt, while a satellite 
trapped in the Kozai resonance has an argument of pericen- 
ter that librates around ^tt or |tt (Kozai 1962; Carruba et al. 
2002). We est imat e the incompleteness in our survey due to 
such orbits in 



Because the equations of motion are symmetric around the 
( = plane, we may further restrict the inclination to the range 
0° < I < 90°. As in the two-dimensional case, we define "pro- 
grade" and "retrograde" in the rotating frame unless otherwise 
noted. Thus retrograde orbits have £ < and prograde orbits 
have £ > at this surface of section 77 = £ = £ = 0, 77 > 0. 

2.2. Numerical orbit integrations 

Even in the two-dimensional case, we expect that the sta- 
ble regions for distant satellites of the giant planets will be 
somewhat different from those derived by Henon (1970) and 
shown in Fig.[U since (i) Henon 's results are based on Hill's 
approximation /i — > 0, while the giant planets have /i in the 
range 0.00096 (Jupiter) to 0.000044 (Uranus); (ii) Henon's 
results assume that the planet orbit is circular, while the gi- 
ant planets have eccentricities between 0.0086 and 0.056; (iii) 
both the satellites and their host planets are subject to pertur- 
bations from the other planets. We must carry out long-term 
numerical integrations of the satellite orbits to assess the in- 
fluence of these effects on the stability region shown in Fig. 

m 

We sample the initial conditions using a fine grid on the 
Henon diagram, with dT = 0.1 and = 0.06. This is shown 
as the dotted grid in Fig. [2] We then convert them to the non- 
rotating (xyz) planetocentric coordinate system where we do 
the integrations of satellite orbits. We require that in the ro- 
tating frame the Sun is always located at the -£ axis, and the 
angular velocity of the rotating frame equals the instantaneous 
angular velocity of the Sun relative to the planet in the non- 
rotating planetocentric frame; thus the angular speed of the 
rotating frame is time- varying if the planet's orbit is eccen- 
tric, and the direction of the ( axis may vary if the planet's 
orbit is perturbed by other planets. We use a unit of length 
\i l l 3 a p and unit of time n~ l to scale the coordinates/velocities 
between the two frames, where a p is taken to be the initial 
semi-major axis of the planet. 

The system to be numerically integrated is composed of 
the four outer giant planets (or sometimes just one of them), 
the Sun, and a satellite around one of the planets; the satel- 
lite is treated as a massless test particle. We use a second- 
order Wisdom-Holman symplectic scheme (Wisdom & Hol- 
man 1991), as implemented in the Swift package (Levison & 
Duncan 1994). Following Nesvorny et al. (2003), we have 
modified the Swift code such that the integration of the plan- 
ets is done in the Jacobi coordinate system while that of the 
satellites is done in the non-rotating planetocentric coordi- 
nate system. We tried different timesteps to optimize between 
speed and accuracy, and found dt = 20 days is short enough 
to produce the correct results with reasonable computational 
cost, for all four planets. 

One potential concern is that the Wisdom-Holman sym- 
plectic scheme, as we have implemented it, is designed for 
nearly Keplerian orbits relative to the planet and might break 
down at large distances from the planet, where the orbits are 
nearly Keplerian relative to the Sun. However, the character- 
istic orbital period at large distances is equal to the planetary 
orbital period, and this is much longer than the orbital pe- 
riods of satellites inside the Hill radius that the integrator is 
designed to follow, so even a crude integrator should work 
well. Moreover, our ability to reproduce the Henon diagram 
(compare Fig. [T] and the lower right panel of Fig. [2]), the long- 
term stability of many of our orbits, and the similarity of th e 
characteristic orbit shapes to those found by Henon (see 33.11) , 
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Fig. 2. — Henon diagrams for Jupiter under various circumstances (see 33.1l for details). The satellite orbital plane initially coincides with Jupiter's. In each 
panel, the dotted region is the grid of initial conditions in the (T,£) plane. The two blank regions in the upper-right corner are forbidden. The filled circles in 
different colors represent the initial conditions of orbits that survive for various times. 



all indicate that even at the largest distances probed here, the 
symplectic integrator seems to work pretty well. As a fur- 
ther check, we have used the Bulirsch-Stoer integrator to fol- 
low satellite orbits around Uranus for 10 6 yr and found almost 
identical results to the Wisdom-Holman integrator. 

We terminate the integration if the distance of the satel- 
lite from the planet exceeds the semi-major axis of the planet 
since at this point the satellite has escaped according to our 
definition at the end of §Q] or if the distance is less than 
the semi-major axis of the outermost regular satellite of each 
planet (being Callisto, Iapetus, Oberon and Triton respec- 
tively), since at this point the satellite lifetime against ejection 
or collision with the regular satellite or the planet is likely to 
be short. Any test particles that cross either of these two radii 
are considered lost. We have experimented with including the 
quadrupole moment J2 of the planet (including the contribu- 
tion from the inner regular satellites) but this has no detectable 



effect on our results. 



3. RESULTS 



3.1. Two-dimensional Henon diagrams 
We first study cases in which the initial velocity vectors of 
satellites lie in the planet orbital plane, i.e., ( = ( = 0. As 
we have described, this is different from Henon 's problem be- 
cause: (i) we do not use Hill's approximation; (ii) planets such 
as Jupiter have non-zero eccentricity; and (iii) there are grav- 
itational perturbations from other planets. 

As an illustration we show how the stable region changes 
under various conditions in Fig. [2 for satellites around Jupiter 
and an integration time 10 6 yr. We consider four situations: 
(a) Jupiter moves on its actual (slightly eccentric) orbit, in- 
cluding perturbations from the other three giant planets (upper 
left); (b) the planar restricted three-body problem, in which 
Jupiter travels on an orbit with its current eccentricity of 0.048 
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Fig. 3. — Two-dimensional Henon diagrams for the four planets. Each orbit is integrated up to 10 8 yr under the gravitational influence of its host planet, the 
Sun, and the other three giant planets. Notations are the same as Fig. [2] 



and the other planets are absent (upper right); (c) the planar 
circular restricted three-body problem, in which Jupiter trav- 
els on a circular orbit with its current semi-major axis (bottom 
left); (d) same as (c) except that the planet mass is 1/100 of 
the Jupiter mass (bottom right). 

By comparing Figs. Q] and [2] it is clear that case (d) in the 
bottom right panel best reproduces the original Henon dia- 
gram; this is not surprising since n ~ 10~ 5 is smallest so 
Hill's approximation is satisfied best, and the other conditions 
assumed in Henon 's problem (circular planet orbit, no other 
planets) are also satisfied. When using the actual Jupiter mass 
in (c), the outer retrograde stable region (i.e., the lower-left 
branch) shifts and shrinks. The overall stability region shrinks 
further — but does not vanish — when Jupiter's orbit is eccen- 
tric as in case (b), and for the most realistic case (a). 

Note that 10 6 yr is only a small fraction of the lifetime of 
the solar system, and possible erosion of the stable region over 
longer times is somewhat indicated by the presence of a few 



red and blue dots in the upper panels of Fig. O indicating 
orbits that are unstable on timescales of 10 4 and 10 5 yr. 

These illustrative calculations show that some Jovian satel- 
lites orbiting well outside the Hill radius can survive for at 
least 10 6 yr, although the stable region is substantially smaller 
than in Hill's approximation to the circular restricted three- 
body problem and appears to erode slowly with time. They 
also show that the stable region is larger (relative to the Hill ra- 
dius) if the planet mass /i is smaller, suggesting that the stable 
regions of the other giant planets may be larger than Jupiter's. 

We now extend these calculations in the following ways: 
(i) we examine satellite orbits around all four giant planets, 
using the actual planetary orbits including perturbations from 
the three other planets; (ii) since stable orbits are found at the 
most negative Jacobi constant (T = -6) examined in Fig. [2 we 
extend the grid of initial conditions to T = -16; (iii) we extend 
the integration time from 10 6 yr to 10 8 yr. 

The results are shown in Fig. [3] There are large regions 
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integrations. In contrast to the results in Fig. [3] Saturn can host stable outer retrograde orbits, and most of the satellites that survive for 10 6 yr also survive for 
10 8 yr around both planets. 



of inner retrograde/prograde orbits that are stable for 10 8 yr. 
For Jupiter, there are a few outer retrograde orbits that survive 
for 10 8 yr; Wiegert et al. (2000) found no orbits that survived 
for > 10 7 yr, but this may reflect their less complete cover- 
age of phase space. For Saturn, the outer retrograde stable 
region completely disappears in less than 10 6 yr, a conclusion 
already reached by Wiegert et al. For Uranus and Neptune, 
in contrast, there is a large stable region of outer retrograde 
orbits remaining after 10 8 yr. We expect that the stable re- 
gions around Jupiter, Uranus, and Neptune will shrink some- 
what further between 10 8 yr and 5 x 10 9 yr, the approximate 
age of the solar system, so we integrated some outer retro- 
grade satellite orbits around these three planets for 10 9 yr. We 
found that about a third of the Jovian orbits and over half of 
the outer retrograde orbits for Uranus and Neptune shown in 
Fig. [3] still survive. Thus it is very likely that Uranus and 
Neptune could still host primordial satellites on such orbits to 
the present time. It is likely, but not certain, that similar satel- 
lites could survive around Jupiter, at least in small volumes of 
phase space. 

The shrinkage of the stable region of the outer retrograde 
branch between 10 6 and 10 8 yr, as well as the lack of stable 
outer retrograde orbits around Saturn, appear to be mainly due 
to perturbations from the other planets. To demonstrate this, 
we ran two 10 8 yr integrations for Saturn only and Uranus 
only. The results are shown in Fig.Hl in this case, Saturn can 
host stable outer retrograde satellites for at least 10 8 yr, and 
there is almost no difference in the size of the stable region 
between 10 6 and 10 8 yr for either planet. 

The stable region around Uranus is larger than the one 
around Saturn in Fig.Hl and the stable regions around Uranus 
and Neptune are larger than the one around Jupiter in Fig. [3] 
These differences are probably caused mostly by their differ- 
ent planet-to-Sun mass ratios /i. As fi increases, the outer 
retrograde stability branch in the lower left of the Henon di- 
agram shrinks, and shifts upward (see Henon 1965; 1970, or 
compare the two lower panels of Fig. [2]). 

We also notice in Figs. [2]l and in the right panel of |4] that 
there is a little tail or branch to the stable region around 



= (-5,-2.5). We suspect this comes from Henon 's pe- 
riodic family #3, which bifurcates from the periodic retro- 
grade orbits at (T, £) = (-2,-1.2) and passes close to the point 
(r,0 = (-5,-2.5) (see Henon 1970, Fig. 13). 

What do these stable orbits look like? In Hill's approxi- 
mation, the stable (outer and inner) retrograde and inner pro- 
grade orbits are generated from the periodic / and g families 
respectively using the terminology of Henon (1969). We show 
examples of stable orbits (i.e., those that survived for 10 8 yr) 
around Uranus in Fig. [5] For each example orbit, we plot the 
instantaneous locations for the first one million years as dots 
with the two stars marking the starting and ending locations. 
We also plot the trajectory for several revolutions. Each or- 
bit is plotted in both the non-rotating planetocentric x-y plane 
(left column) and the rotating £-77 plane (right column). In the 
rotating frame, the inner prograde orbit (top panel) is elon- 
gated along the Sun-planet axis while the inner retrograde 
orbit (middle panel) is elongated perpendicular to the Sun- 
planet axis. The outer retrograde orbit (bottom panel) is also 
elongated perpendicular to the Sun-planet axis and oscillates 
about the planet as an ellipse with an axis ratio of approxi- 
mately 2:1 (compare Fig. 11 of Henon 1970), as one would 
expect from epicycle theory. 

Note that the stable retrograde orbits with £ close to 
-0.6934 in Fig. [3] regularly cross the Hill radius. Hence such 
orbits are missed by the surveys of both Wiegert et al. and 
Nesvorny et al., who terminate their integrations if r < ru or 
r > rn, respectively. 

3.2. Three-dimensional Henon diagrams 

We now extend the initial conditions in 33.1l to three dimen- 
sions by including the scaled vertical coordinate f. As dis- 
cussed in 32.11 we consider a surface of section r\ = ( = £ = 0, 
77 > at t = 0. Similar to the two-dimensional case, we 
sample the initial conditions using a fine grid in the (T,Q 
plane, and use equation © to generate initial velocities. We 
choose a sequence of inclinations in the rotating frame, / = 
15°,30 o ,45 o ,60 o ,75°. Each satellite orbit is then integrated 
for 10 8 yr along with the four giant planets and the Sun. 
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Fig. 5. — Examples of stable orbits around Uranus. All orbits initially lie in the orbital plane of Uranus. Dots are instantaneous locations for the first 10 6 
yr, plotted at intervals of 100 yr, with the green and red stars marking the starting and ending locations. We also plot a few revolutions as red curves. The left 
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Hill sphere. Upper: an inner prograde orbit; middle: an inner retrograde orbit; bottom: an outer retrograde orbit. 
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The general behavior when incorporating inclination is the 
erosion of stable regions in the Henon diagram. As an ex- 
ample, we show the results for Uranus in Fig. [6] The outer 
retrograde orbits quickly become unstable when the initial in- 
clination exceeds ~ 20° . The inner retrograde stable region 
erodes with increasing inclination and disappears at / > 75°. 
The inner prograde stable region can survive even at / w 75°. 
This asymmetry between inner retrograde and prograde or- 
bits is due in part to the definition of inclination in the rotat- 



ing frame 3 . However, even when inclination is defined in the 
non-rotating frame such an asymmetry may still be present 
(see Cuk & Burns 2004 for a discussion of the dynamical rea- 
sons for the asymmetry). 

To separate the destabilizing effects of inclination from the 
effects of perturbations from other planets, we also ran these 

3 When translated into the non-rotating planetocentric frame, the inclina- 
tions of "prograde" ("retrograde") orbits are actually smaller (larger) than 
in the rotating frame. As we already noted in Fig. [T] under certain circum- 
stances retrograde orbits in the rotating frame can even be prograde in the 
non-rotating frame. 
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three-dimensional simulations for Uranus without the other 
planets. We found that for I < 20°, perturbations from other 
planets do play a major role in eroding the region of stable 
outer retrograde orbits, as illustrated by comparing the upper 
left and lower right panels of Fig. [6l which show the Henon 
diagram for 7=15° with and without the other planets. How- 
ever, for I = 30° , 45° , 60° , 75° the results for Uranus alone are 
almost identical to the realistic case which includes perturba- 
tions from the three other planets. This result suggests that 
bound retrograde orbits outside the Hill sphere may not exist 
at all for I > 20° . This is expected because the Coriolis force, 
which stabilizes outer retrograde orbits by bending their tra- 
jectory towards the planet in the rotating frame, is reduced 
when the inclination angle I increases. This shrinkage of the 
stable outer retrograde branch with inclination is already no- 
ticed by comparing the right panel of Fig. [4] (I = 0) and the 
bottom-right panel of Fig. 0(7 =15°). 

Note that the inclinations I are planetocentric and measured 
in the rotating frame. For low-inclination orbits along the 
outer retrograde branch they can be converted to heliocentric 
inclinations i using the approximate formula for small e 



i « 2el , 



(10) 



where e is the heliocentric eccentricity. Thus our stability 
region I < 20° corresponds roughly to i < 4° for e w 0.1, 
the typical heliocentric eccentricity of surviving satellites 
for Uranus and Neptune (see Fig. [12]); this is in reasonable 
agreement with Wiegert et al.'s estimate that most of their 
long-lived orbits had i < 2°, especially considering that our 
sampling of initial conditions is more complete than theirs. 
Mikkola & Innanen (1997) and Mikkola et al. (2006) esti- 
mate analytically that the outer retrograde orbits are unstable 
if i > e (for a circular planet orbit), which implies instability 
if I > 30°. Our own results show that all the test particles with 
initial positions outside the Hill sphere cross the escape radius 
a p well before 10 4 yr for I > 30°. 

As described at the end of 32.11 a limitation of these results 
is that the Henon diagram we have used will not display or- 
bits trapped in a Kozai resonance, or other stable orbits whose 
argument of pericenter does not periodically pass through 
or 7r. To estimate the contribution of such orbits, we have 
constructed a different set of three-dimensional Henon dia- 
grams in which the initial conditions are changed from our 
usual choice rj = £ = ( = 0,7] > to = £ = £ = 0,7j>0 (i.e., 



when the orbit is at its maximum height above the planet's 
orbital plane, rather than crossing the orbital plane). In this 
case, we define the initial inclination 7* in the rotating frame 
by 

tan/* = j : . (11) 

£ f=0 

The results are shown in Figure [71 which should be compared 
to Figure [6) Each point in either set of Henon diagrams corre- 
sponds to a unique orbit, but orbits appearing in the Henon 
diagrams of one figure at a given value of (T,^,7) may or 
may not appear in the other figure, where they will have the 
same value of T but possibly different values of £ and incli- 
nation. The stable regions are somewhat larger in Figure [7] 
at a given inclination — for example, a few outer retrograde 
satellites survive for 10 8 yr at 7* = 30° — but the conclusions 
described above are not significantly altered. In the follow- 
ing discussion, we neglect stable orbits that do not appear in 
our fiducial Henon diagrams (i.e., using the surface of section 
77 = £ = C = 0,77>0); thus we may slightly underestimate the 
size of the stable regions. More discussion on orbits trapped 
in the Kozai resonance inside the Hill sphere can be found in 
Carruba et al. (2002). 

3.3. Spatial stability regions 

We now project the phase-space volume that hosts stable 
orbits onto coordinate space, to explore where stable satellites 
might be found. 

We plot the positions of stable orbits in the two-dimensional 
plane with coordinates [(x 2 + v 2 ) 1 / 2 ,z] (non-rotating frame) or 

[(£ 2 + ?? 2 ) 1 / 2 , C] (rotating frame). Prograde and retrograde or- 
bits are plotted separately on the left and right sides of a given 
figure panel, where "prograde" and "retrograde" are defined 
in the frame used. We plot the position of each stable (up to 
10 8 yr) point in the Henon diagram at uniformly spaced times 
(every Myr) between 5 x 10 7 and 10 8 yr in Figs. [8] (Jupiter) to 
GJ] (Neptune). 

The stability regions within the Hill sphere are very similar 
to those shown in Figs. 9-12 of Nesvorny et al. (2003), though 
slightly larger because we show instantaneous position rather 
than semi-major axis. For Jupiter and Saturn, the stable pro- 
grade orbits generally extend to ~ 0.5rn; the stable retrograde 
orbits can extend further to ~ 0.7rn, and nearly coplanar ret- 
rograde orbits even extend to ~ rn for Jupiter. For Uranus and 
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Neptune, both the prograde and retrograde stable orbits can 
extend a little bit further relative to the Hill sphere. No sta- 
ble orbits exist at high latitudes, presumably because of Kozai 
oscillations (Kozai 1962; Carruba et al. 2002; Nesvorny et al. 
2003). 

It is also notable that there are stable regions beyond the 
Hill sphere for Jupiter, Uranus and Neptune, as we dis- 
cussed in previous sections. This is particularly the case for 
Uranus and Neptune. Most of these distant stable satellites 
are concentrated close to the orbital plane of the planet as 
for Jupiter, although the appearance of a very thin layer in 
Figs. [8] (Jupiter)- Qj] (Neptune) is somewhat an artifact of 
the coarse sampling of inclinations in our initial conditions 
(/ = 15°,30 o ,45 o ,60 o ,75°). Stable satellites can be found as 
far as ~ 5rn from Jupiter and even ~ lOrn for Uranus and 
Neptune, and as high as 2.5r H above the orbital plane for the 
latter two planets. In the following section we discuss briefly 
the strategy of searches for these distant satellites. 

4. DISCUSSION AND CONCLUSIONS 

We have conducted a systematic survey of the stable regions 
of satellites around giant planets in the solar system, using nu- 
merical orbital integrations that include gravitational pertur- 
bations from the other planets. We confirm previous results 
for satellites within the Hill sphere: stable retrograde satel- 
lites can exist further out than prograde satellites (e.g., Henon 
1970; Hamilton & Krivov 1997; Nesvorny et al. 2003); and 
stable orbits cannot exist at high inclinations (e.g., Carruba et 
al. 2002; Nesvorny et al. 2003). 

We also confirm and extend the conclusions of Wiegert et 
al. (2000) that distant retrograde satellites ("retrograde" as de- 
fined in the rotating frame) can survive well beyond the Hill 
sphere for at least 10 8 -10 9 yr, and probably for the lifetime of 
the solar system. Uranus and Neptune are the most promising 
host planets for such distant satellites, since their stability re- 
gions have the largest extent (e.g., Fig. [TQlfTTb . Jupiter has a 
smaller stable region (Fig. [5]), and Saturn appears to have no 
stable regions beyond the Hill radius (Fig. [9]). 

Remarkably, there is a gap between the inner and outer sta- 
bility zones for retrograde satellites, extending from about rn 
to 2rn, in which almost no stable orbits exist. 

To check whether any of the proposed distant satellites have 
already been discovered as Centaurs, we take the positions 

4 http://cfa-www.harvard.edu/iau/mpc.html 



and velocities of the known Centaurs from the IAU Minor 
Planet Center 4 that have planetocentric distances smaller than 
the semi-major axes of each of the four giant planets at the 
last observed epoch. There are 31 Centaurs (Jupiter 1; Uranus 
16; Neptune 14) that satisfy the criterion. We numerically 
integrated these objects along with the Sun and giant planets 
for 10 8 yr, but none of them survived as satellites according 
to the definition in §Q] 

Searches for satellites far beyond the Hill radius can be 
carried out either with dedicated deep, wide-angle surveys 
around the giant planets, or through all-sky surveys such as 
Pan-STARRS and LSST. The most promising search areas 
are close to the orbital plane of the planet, since only low- 
inclination orbits survive (e.g., compare Fig. [3] and Fig. [6]). 
In Fig. [T21 (left panel) we show the heliocentric angular ve- 
locity of the stable outer retrograde satellites as a function of 
angular distance from the planet as viewed from the Sun, for 
Jupiter (black), Uranus (blue) and Neptune (red) respectively, 
sampled every Myr between 5 x 10 7 and 10 8 yr. 

In the middle and right panels of Fig. [321 we show his- 
tograms of the difference in heliocentric semi-major axis from 
their host planet and heliocentric eccentricities for the stable 
outer retrograde satellites. These distributions can be used to 
cull a large sample for potential satellites. Once a candidate 
is identified with reliable orbit elements, a long-term orbital 
integration should be run to confirm its satellite nature. 

The discovery and characterization of satellites beyond the 
Hill sphere would provide rich information about the early 
formation of the solar system. Fabry cky (2008; also see Ko- 
rtenkamp 2005) recently performed simulations of capture of 
neighboring planetesimals from the circumstellar disk during 
slow planet growth, and found that such distant satellites are 
a natural outcome for Uranus and Neptune. Thus an inven- 
tory of this potential population of bodies would enhance our 
understanding of the formation of planets and their satellites 
in the early solar system, and the properties of the primordial 
planetesimal disk. 
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